#! /bin/sh
# This input file reproduces the output in Juhlin (1995) Figure 8b.
# Note that the source in Table 1 in Juhlin (1995) is given as g(t), but should be f(t)
# be f(t), i.e. the derivative of the Gaussian (dg)
# Since the dg source waveform is defined by its avg frequency content
# favg should be set to 17.077 Hz to correpsond to alpha=5000 in
# Juhlin (1995). Also the time step, dt, is not 0.8 ms, but 0.5 ms
#
# The unfiltered output contains high frequency numerical noise (I do not
# quite understand why) which is nicely removed usinf a highcut filter
#
# Note that the geometry of the model in Figure 8a is very schematic 
# and does not correspond to the model actually used in producing the 
# synthetics (the things we can get away with in our publications).
#
#suea2df dt=0.0005 lt=1.2 zmax=3500 xmin=-00 xmax=4000 dx=5 verbose=1 \
#sx=0 sz=20 favg=17.077 ts=0.08 wtype=dg \
#snfile= hsz=2000 vsx=2000 snaptime=0.19,0.59,0.99 efile=mod \
#bc=2,10,10,10 qsw=0 asw=1 msw=1 \
#xl=0,3500 \
#0 \
# zl=0,5,500,505 \
# vpl=3320,2898,2898,3320 \
# vsl=1398,1398,1398,1398 \
# rhol=2250,2250,2250,2250 \
# epl=0,0.2566,0.2566,0\
# dsl=0,0.4025,0.4025,0\
# anl=0,-30,-30,0 \
#3500 \
# zl=2000,2005,2500,2505 \
# vpl=3320,2898,2898,3320 \
# vsl=1398,1398,1398,1398 \
# rhol=2250,2250,2250,2250 \
# epl=0,0.2566,0.2566,0\
# dsl=0,0.4025,0.4025,0\
# anl=0,30,30,0 \
#>out2

modelfile=model3.unif           # input model file for unif2aniso
ninf=3                          # number of interfaces (surface counts)
x0=0,0,0                        # x-position(s) for  vp00,vs00,rho00, etc.
z0=0,5,500,505 			# z-position(s) for  vp00,vs00,rho00, etc.
nz=700                          # size of z (depth) dimension  of model
nx=800                          # size of x (horizontal) dimension of model
dz=5                          # increment in z direction
dx=5                          # increment in x direction
vp00=3320,2898,2898,3320        # P-wavespeed(s) at (z0,x0)
vs00=1398,1398,1398,1398        # S-wavespeed(s) at (z0,x0)
rho00=2250,2250,2250,2250       # density(s) at (z0,x0)
eps00=0,0.2566,0.2566,0	        # epsilons(s)
delta00=0,0.4025,0.4025,0	# delta(s) at (z0,x0)
gamma=0,30,30,0 	        # gammas(s) at (z0,x0)
method=linear			# linear, mono, akima, spline  interpolation
                                #        method in unif2aniso
dt=0.0005                       # time sampling interval
lt=1.2	                        # latest time modeled
fx=0                            # first x value
verbose=1                       # =1 chatty, =0 silent
snfile="snaps.su"               # output file for snapshots
rhofile="rho_file"              # input file of densities, if not set rho=1000
hsz=0		                # z-position of horizontal line of geophones
vsx=20		                # x-position of vertical line of geophones
snaptime=.04,.08,.06,.12,.14,.16,.18,.2,.22,.24,.26,.28,.30,.32,.34,.36,.38,.40,.42,.44,.46,.48,.50,.52,.54,.58,.60,.62,.60,.66,.68,.70,.72,.74,.76,.78,.80,.82,.84,.86,.88,.90.,92,.94,.96,.98,1.00,1.02,1.04,1.06,1.08,1.10
# times of snapshots
bc=2,10,10,10                   # boundary conditions
qsw=0                           # =1 put in attenuation
asw=0                           # =1 anisotropy
sx=10				# x-position of sources
sz=20				# z-position of sources
favg=24.15                      # average frequency
ts=.05                          # source duration
wtype=dg                        # waveform type


# build stiffness and density files
unif2aniso < $modelfile ninf=$ninf x0=$x0 z0=$z0 nz=$nz nx=$nx \
dx=$dx dz=$dz vp00=$vp00 vs00=$vs00 rho00=$rho00

xbox=10
ybox=10
nxplot=`bc -l <<-END
	scale=0
	$nx / 3 
END`
nzplot=`bc -l <<-END
	scale=0
	$nz / 3 
END`

# the files c11_file c13_file c15_file c33_file c35_file c55_file rho_file
# are generated by unif2aniso
# transpose stiffness and density
for i in c11 c13 c15 c33 c35 c55 rho
do
	echo $xbox $ybox

	ximage <  ${i}_file n1=$nz n2=$nx perc=99 xbox=$xbox \
		ybox=$ybox  wbox=$nxplot hbox=$nzplot  legend=1 \
		title=" ${i} parameter file  "	&

	xbox=`expr $xbox + 110 `
	ybox=`expr $ybox + 5 `

        mv ${i}_file tmp.file
        transp n1=$nz < tmp.file > ${i}_file
done

rm tmp.file


# run suea2df
suea2df dt=$dt lt=$lt nz=$nz fx=$fx nx=$nx dx=$dx dz=$dz verbose=1 \
snfile=$snfile rhofile=$rhofile hsz=$hsz vsx=$vsx snaptime=$snaptime \
bc=$bc qsw=$qsw asw=$asw sx=$sx sz=$sz favg=$favg ts=$ts wtype=$wtype \
>out2


n2=`bc -l <<-END
      scale=1
       $nx * 2
END`

echo $n2


suxmovie < snaps.su n1=$nz n2=$n2 clip=1e-13 loop=1 title="snapshots horizontal vertical " width=$nxplot height=$nzplot  sleep=70000 &



# shot gathers from a horizontal line of geophones
suximage <  hs.su xbox=0 ybox=400 wbox=$nxplot hbox=$nzplot  perc=99 title="  shot gathers  vertical horizontal " &
suximage <  vsp.su xbox=400 ybox=400  wbox=$nxplot hbox=$nzplot  perc=99 title=" VSP  vertical horizontal " &

exit 0
